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Abstract. We present an algorithm for the computation of unbiased Green functions and self¬ 
energies for quantum lattice models, free from systematic errors and valid in the thermodynamic 
limit. The method combines direct lattice simulations using the Blankenbecler-Scalapino-Sugar 
quantum Monte Carlo (BSS-QMC) approach with controlled multigrid extrapolation techniques. 
We show that the half-filled Hubbard model is insulating at low temperatures even in the weak- 
coupling regime; the previously claimed Mott transition at intermediate coupling does not exist. 


1. Introduction 

Numerous studies of the Hubbard model [l have greatly enhanced our understanding of strongly 
correlated electron systems within the last four decades. Specialized methods, namely the semi- 
analytic Bethe ansatz and the density matrix renormalization group (DMRG) |2|, yield reliable 
high-precision results (only) in the case of one spatial dimension. Conversely, the dynamical 
mean-held theory (DMFT) 3, [4] provides deep insight in the limit of high dimensionality. 
However, the intermediate regime of two (and three) dimensions is less well understood, e.g., 
with respect to pseudogap physics |5j|7| and high-T c superconductivity. In this paper, we will 
address another open question: the nature of the Mott metal-insulator transition (MIT) |8j|9] of 
the half-filled Hubbard model in two dimensions (D = 2). 

As indicated by shaded regions in Fig. [lj the Hubbard model is insulating at half filling (1 
electron per site) and sufficiently strong coupling, in any dimension, while it is metallic at weak 
coupling - as long as the translational symmetry is not broken, e.g., by magnetic order. For 
D = 3, one finds antiferromagnetic (AF) ordering |10|, which implies insulating behavior, below 
the Neel temperature TV (gray dash-dotted line). Such long-range order is ruled out by the 
Mermin-Wagner theorem in D = 2 (at temperature T > 0). Single-site DMFT predicts a MIT 
with critical interactions 2.3 < U c < 2.9 (pink dotted line). This transition line is shifted to 
weaker interactions within cluster DMFT (CDMFT), which includes short-range correlations 
on small clusters (violet solid line) in. However, the low-T metallic phase implied by the 
CDMFT result is incompatible with several earlier studies highlighting effects of “short-range 
antiferromagnetism” [7,12, 13 . To clarify the situation, we apply the Blankenbecler-Scalapino- 
Sugar quantum Monte Carlo (BSS-QMC) [l4] approach with controlled multigrid extrapolation 
techniques to selected points in the phase diagram (green crosses). We will show that these points 
are separated by a MIT line (or narrow crossover) and conclude that the 2D Mott physics is 
quite similar to the 3 D case [15], in spite of the Mermin-Wagner theorem. 
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Figure 1. Schematic representation of the 
MIT in the Hubbard model at half filling on 
a square lattice. Green crosses: parameters 
selected for our unbiased QMC study. Gray 
dash-dotted line: Neel temperature for the 
cubic lattice (rescaled). 



Figure 2. Extrapolation of BSS-QMC 
Green functions to At — > 0. Analytic con¬ 
tinuation (MEM, black line) regularizes the 
results of pointwise extrapolations (symbols) 
and yields stable results (gray line) even in 
the case of noisy raw data (thin color lines). 


2. Model and Methods 

The single-band Hubbard Hamiltonian is given by 

h=-1 + u Yl (%- (W - > (!) 

(vV i V J V 7 


with number operators h la = c| c- , next-neighbor hopping t, and local interactions U. In this 
paper, we set t = 0.25 and consider finite-size clusters with periodic boundary conditions. 

The BSS-QMC algorithm is based on a Trotter-Suzuki decomposition of the partition function 


Z = Tr ( e ~^ Ht+Hu « Z At = Tr e _ATfft e _Ar/fu ^ , (2) 

where Hu ( Ht ) corresponds to the interaction (kinetic) term in ([!]), /3 = 1/T is the inverse 
temperature (kb = 1) and A denotes the number of time slices (A = /3/A r). A discrete Hubbard- 
Stratonovich transformation simplifies the interaction term, which is quartic in the fermionic 
operators, to a quadratic form and a coupling to an auxiliary Ising field h. As usual for quadratic 
Hamiltonians, the fermionic degrees of freedom can be integrated out and the partition function 
is expressed by determinants of matrices: 
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The sum in ([3]), and finally the computation of all observables, is performed by Monte Carlo 
methods. The main results of the BSS-QMC algorithm are Green functions Gij(r{) on a discrete 
imaginary-time grid {r; e [0, (3]} for each pair of real-space coordinates {i,j} on the given lattice. 
The involved matrix operations for the calculation of the determinants in ([3]) lead to a scaling 
of O(ANc), which restricts the algorithm to relatively small clusters ( N c < 400). 

To arrive at reliable conclusions regarding the Mott transition using BSS-QMC, one has first 
to eliminate all systematic errors (finite-size as well as Trotter errors) and, secondly, calculate 
the self-energy on the Matsubara axis, E(*w n ), from the imaginary-time Green function in a 
stable manner. These steps are described in the following. 



















Figure 3. (a) Influence of statistical noise in the raw data G(t) on the self-energy Y,[iu n ) 
and stabilization by MEM regularization. Thin lines: noisy data, thick lines and symbols: high 
precision results, (b) Impact of Trotter discretization At on the Matsubara self-energy. 


16-19 , we will sketch the scheme only 


3. Elimination of systematic errors 

As the elimination of Trotter errors is well established 
briefly and focus on the problem specific and physically more relevant finite-size analysis and 
present a stable method to calculate unbiased Green functions, even if the raw data are noisy. 

Trotter bias - The raw BSS-QMC Green functions (thick colored lines in Fig. [2]) do not 
only live on different imaginary-time grids {r*}, depending on the chosen discretization At, but 
are also shifted with respect to each other (and to the exact solution). After aligning them on 
a common fine grid [201 the Trotter errors can be eliminated with high precision by piecewise 
extrapolations of At —> 0 (circles) 1191. However, some fluctuations remain inevitably, especially 
when using noisy raw data (using a factor of 10 less QMC sweeps; thin colored lines). These can 
be greatly reduced by regularization via a maximum entropy method (MEM), which computes 
corresponding spectral functions via the inversion of 


/ 


G(t) = — du A(u) 


1 + e~“P 


( 4 ) 


Note that the intermediate spectra A(u) are used only for producing continuous and smooth 
Green functions G(t ) via Q; in this combination, the procedure is stable. After this step, 
the results based on high- or low-precision data (thick black line vs. gray line) are hardly 
distinguishable and can both be used for stable computations of self-energies. 

Self-energy - These quasi-continuous G(r) can be reliably Fourier transformed to the 
Matsubara axis: 

r? 

G(iu n ) = dr G(t) e . (5) 

Jo 

A quantity of great interest for the analysis of the Mott transition is the imaginary part of the 
momentum-resolved self-energy, which is connected to the BSS-QMC Green function G and the 
non-interacting Green function Q via a Dyson equation [21 

E k (Jojji) — (ic^TT,) (> k ( ioj n ) — ioj n T /r e k G k , (6) 


where fi is the chemical potential and the dispersion of the noninteracting problem. As 
shown in Fig. [3] (b), the MEM regularizing procedure leads to estimates (crosses and black line) 
of the imaginary part of the self-energy (the real part vanishes at k = (tt, 0) due to particle-hole 
















Figure 4. Finite-size scaling of Matsubara self-energy at U = 0.5, k = (7r, 0). Finite- 
size BSS-QMC data (open symbols and broken colored lines), extrapolated BSS-QMC results 
in the thermodynamic limit (circles and black bold solid line), and Dr A data (gray line) 
versus Matsubara frequency oj n at T = 0.01 (a) and T = 0.025 (c); also shown are 
momentum-independent single-site DMFT results (thin black lines). (b)+(d) Finite-size BSS- 
QMC (symbols) data for the first three Matsubara frequencies versus inverse system size plus 
extrapolations in linear order in L -2 (thin lines) and quadratic order (thick lines). 


symmetry) which are smooth even at higher frequencies. In contrast, the direct results (green 
circles) fluctuate significantly even when based on high-precision raw data (and even more so for 
noisy raw data; thin dashed green lines). Note that both methods agree very well at the lowest 
Matsubara frequencies, where the QMC predictions are most reliable. 

In Fig. [3] (b) the influence of the Trotter discretization is shown. The Trotter biased self¬ 
energies (broken colored lines) deviate from the exact result (black solid line) for large values of 
Ar, in particular at the lowest frequencies. However, in the parameter regime explored in this 
study, we can afford small enough values of Ar so that the elimination of Trotter errors is less 
essential than the finite-size extrapolation to be discussed below. 

Finite-size extrapolation of the self-energy - Using the methods described above, unbiased 
and stable estimates of the Matsubara self-energy have been obtained for square L x L clusters 
with sizes ranging from 8 x 8 to 16 x 16 (with a computational effort that varies as = T 6 , i.e., 
by a factor of 2 6 = 64 in this range). Results at the “anti-nodal” |22| momentum point k = (r, 0) 
[throughout the paper, unit lattice spacing a = 1 is assumed] are shown (colored symbols and 
broken lines) in Fig. |4](a) for the low temperature T = 0.01 and in Fig. |4](c) for the elevated 
temperature T = 0.025, respectively. Evidently, the finite-size effects are enormous: while the 
smallest systems (8x 8, triangles) show insulating behavior, i.e., a strong enhancement of T,(iu; n ) 
towards small frequencies, at both temperatures, this tendency is reduced with increasing L at 
T = 0.01 [Fig. |4](a)] and completely eliminated at T = 0.025 [Fig. |4](c)] . Obviously, careful 
extrapolations are needed for reliable predictions in the thermodynamic limit. 

These are, indeed, possible, as shown in Fig. [4|(b) for T = 0.01 and in Fig. [4](d) for 
T = 0.025, respectively, for the three lowest Matsubara frequencies (where the finite-size effects 
are largest): as a function of L~ 2 , the finite-size results (symbols) can be fitted with second-order 
polynomials (thick lines) with high precision; as these fit functions have small curvatures, linear 
extrapolations (thin lines) to the thermodynamic limit (i.e., L~ 2 —> 0) deviate only slightly from 
the quadratic ones. We use these deviations as error bars and the arithmetic average of both 
extrapolation procedures as final result, as indicated for E(zu;o) by a black symbol in Fig. |4](b) . 
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Figure 5. Real-space spin correlation function, (a) The correlation length for T = 0.025 (red 
squares) from BSS-QMC calculations for a 16 x 16 lattice is estimated by an exponential fit (red 
solid line). For T = 0.01 (green circles) £ is much larger than the extend of the lattice. Right 
panel: The corresponding estimations from DTA data [251 for T = 0.01 (b) and T = 0.025 (c) 
confirm the BSS-QMC results. 


4. Mott physics and correlation lengths 

These final QMC results [black symbols and lines in Fig. Qa) and Fig. Qc)] show that the 
character of the system changes drastically between the selected phase points (crosses in Fig. 
[I]): while the QMC self-energy indicates a metallic phase at T = 0.025, very similar to the 
DMFT solution [thin black line in Fig. [ijc)], the low-temperature phase is clearly insulating 
- and completely unlike the paramagnetic DMFT solution [thin black line in Fig. Qa)]. In 
contrast, the QMC results are strikingly similar to the DTA predictions (gray dash-dotted lines), 
especially at low T. The dynamical vertex approximation (DTA) [23,24j is a diagrammatic 
extension of DMFT which works directly in the thermodynamic limit (without any need for 
finite-size extrapolations); for the full set of DTA results see [25]. QMC and Dr A together show 
convincingly that the suspected phase transition (dotted line in Fig. [I]) is, indeed, present. 

The question remains how exactly this metal-insulator transition (or crossover) is related to 
magnetic properties. After all, long-range AF order is excluded, in D = 2, by the Merrnin- 
Wagner theorem. However, as shown in Fig.[5](a), the spin correlations decay much more slowly 
at T = 0.01 (circles) than at T = 0.025 (squares). More precisely, the correlation length at the 
elevated temperature can be estimated (using 16 X 16 QMC data) to £ « 3.7 while it is clearly 
too large (£ S> 30) at T = 0.01 to allow precise determination using QMC. For this observable, 
DTA predictions, shown in Fig. [5][b) and Fig. [5](c) are more reliable, yielding correlation lengths 
of £ ~ 4 at the higher temperature, in excellent agreement with QMC, and £ > 300 at T = 0.01, 
also consistent with QMC. Using DTA we could also show that the temperature dependence of 
£ changes its character at the Mott transition (see |25|). 

Conclusion and acknowledgements 

In spite of its simplifications, the Hubbard model remains a challenging target for theorists, 
especially in the interesting case of two spatial dimensions. All available methods have certain 
limitations, many of which are of particular concern in the nonperturbative intermediate- 





































coupling regime at low temperatures and in the presence of longer-range correlations. Therefore, 
a single method can hardly yield authoritative results; some predictions may even be completely 
off (such as the CDMFT transition line in Fig. [I]). 

In the study presented in this work (with focus on the QMC methodology; for the full 
story, see |25j), we have applied two methods (an unbiased variant of BSS-QMC as well as 
the dynamical vertex approximation) with completely different characteristics. Only the near¬ 
perfect agreement between both sets of results makes the predictions truly authoritative (and 
validates technical choices on either side). We have settled one important question at half filling, 
namely the character of the Mott metal-insulator transition on the square lattice: it is driven 
by exponentially long-ranged [25] antiferromagnetic correlations, which act similarly to the AF 
long-range order in the cubic case, and is not connected to a quantum-critical point at U > 0. 

The same methodology may be useful in other parameter ranges, e.g., for frustrated or 
anisotropic Hubbard models, doped systems, or multi-band models. It might also be worthwhile 
to compute transport properties, within bubble approximation or beyond, in order to observe 
the metal-insulator transition even more directly (than via the self-energy). 

We acknowledge support from the research unit FOR 1346 of the German Research 
Foundation (DFG) and the graduate school GSC 266. 
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